##############################
# Replication code for Harbers and Ingram  
# PSRM paper, revision, and revision #2
# Previous versions presented at APSA 2016;Southwest Mixed-Methods Workshop 2016; 
# SPSA 2017; APSA 2017; APSA 2018
# Matthew C. Ingram
# University at Albany, SUNY
# contact: mingram@albany.edu
# Last updated: 20181206
##############################

# NOTE:
# Run "harbersingram_psrm_setup.R" script prior to running this main script
# see README file

#################################################################
# 
# READING/IMPORTING AND WRITING/EXPORTING DATA (shapefile, .shp)
#
#################################################################

# If already ran this main script and are returning 
# to data after analysis, just load workspace (.RData file)
# if no data yet, then skip this line
# load("./data/working/harbersingram_smr_20181109.RData")

# If need to get shapefile info use getinfo.shape
# getinfo.shape("NAT.SHP")
# load shapefile; this is original data from Baller et al. (2001)
shp <- readOGR(dsn='./shapefiles', layer="NAT")

# note: sometimes dbf file will not load if extension is .DBF; if this happens, try renaming to lower case .dbf

# load state shapefile; this is from U.S. Census Bureau
# just has state boundaries for contiguous states
shp.st <- readOGR(dsn='./shapefiles', layer="tl_2014_us_state_contig")

# check shp files
#plot(shp)
#plot(shp.st, border="red", add=TRUE)
#summary(shp)

# extract data from shapefile
data<-data.frame(shp)
#names(data)
#head(data)

#####################################################
#
# OLS global model, with regional dummy
#
#####################################################

ols1 <- lm(HR90 ~ RD90 + PS90 + MA90 + DV90 + UE90 + SOUTH, data = data)
summary(ols1)
AIC(ols1)

#####################################################
#
# OLS global model, without regional dummy
#
#####################################################

ols2 <- lm(HR90 ~ RD90 + PS90 + MA90 + DV90 + UE90, data = data)
summary(ols2)
AIC(ols2)

# will need residual from ols2 later; ols2$residuals
# so, write to shapefile that will be imported into python
# and then re-imported in R below

shp$ols2res <- ols2$residuals
writeOGR(shp, dsn='./shapefiles', layer='NAT2', driver="ESRI Shapefile")

###################################################################
#
# Generate grid reference
#
###################################################################

# We need a grid reference for the location of each observation.
# we create this location using the coordinates of the centroids:

# need to merge X-Y coords with shapefile and then re-load shapefile and attach data
#XYgridtable <- cbind(data$XCNTLONG, data$YCNTLAT)
#coords <- XYgridtable
# or simply:
coords <- coordinates(shp)

###################################################################
#
# Generate Spatial Weights
#
###################################################################

# Here, only calculate 10 nearest neighbors (10nnb) 
# this is what is used in running example in paper (Baller et al.)

row.names(shp)
IDs<-row.names(as(shp, "data.frame"))
w_10nnb<-knn2nb(knearneigh(coords, k=10), row.names=IDs)

# make symmetric
wsym_10nnb <- make.sym.nb(w_10nnb)
wsym_10nnb
# and make it a list object
listwsym_10nnb <- nb2listw(wsym_10nnb, style="W")
listwsym_10nnb

# create shortcut to symmetric list object
w <- listwsym_10nnb

# check
#plot(shp)
#plot(w, coords, col="red", add=T)


########################################################################
#####
##### Test global spatial autocorrelation (dependence)
#####
########################################################################

# global dependence of DV
#moran.test(shp$HR90, w)
# global dependence of residuals from ols2
#moran.test(ols2$residuals, w)

# quick local moran plots

#moran.plot(shp$HR90, w)
#moran.plot(ols2$residuals, w)

# local moran plots with ggplot2

###############################################
### Cleaner moran plots with ggplot2
###############################################

###################################
#
# Figure 1
#
###################################

########################################################
# Local Moran and LISA maps with permutation method in Python (pysal)
# Note: permutation analysis done separately in python

# call python script using source_python from reticulate package
source_python("./code/harbersingram_psrm3_20181206_LISAs.py")

# for details, see python script in "code" folder
# results saved to file from python and now imported here
########################################################

# re-create shortcut to symmetric list object in R
# in replication at PSRM, different "w" was being pulled in from python
w <- listwsym_10nnb


# load shapefile with all data from python analysis 
shp2 <- readOGR(dsn='./shapefiles', layer="NAT3-8")
#plot(shp2)

shp <- shp2

shp@data$id = rownames(shp@data)
shp.points <- fortify(shp, region="id")

# if get error about gpclibPermit not being TRUE, try re-installing rgeos:
#install.packages("rgeos")
# if that does not work, try:
#install.packages("gpclib", type="source")
#library(gpclib)
#gpclibPermit()

shp.df = join(shp.points, shp@data, by="id")
#shp.df[shp.df==0] <- NA

#names(shp.df)

###############################
#
# values centered on zero
#
###############################

DV <- shp$HR90
quadrant <- vector(mode="numeric",length=nrow(data))
# center on mean
cDV <- DV - mean(DV) 
lagDV <- lag.listw(w, DV)
# center on mean
clagDV <- lagDV - mean(lagDV)

cl <- as.factor(shp$ls.py.clper)
sig <- shp$ls.py.pperm

dat2 = data.frame(y = cDV, Wy = clagDV, cl=cl, sig=sig, 
                  meany=mean(cDV), meanWy=mean(clagDV))

#############################
# Figure 1: moran plot, b&w 
#############################

tiff(file="./figures/fig1_moranploty_bw_cen.tif", height=6, width=6, units="in", res=300) # could also do pdf, jpeg, bmp, tiff
ggplot(dat2, aes(x = y, y = Wy)) +
  geom_point(colour="black", pch=21, size=3,
             aes(fill = factor(cl))) +
  scale_fill_manual(name = "Cluster",
                    values = c("0" = "white",
                               "1" = "grey0",
                               "2" = "grey30",
                               "3" = "grey60",
                               "4" = "grey90"),
                    labels = c("n.s.", "high-high", "low-low", "low-high", "high-low")) +
  #labs(x="", y="", title="LISA Clusters (DV: HR90)") +
  geom_smooth(method = "lm", se=F, colour="black", size=.7) + 
  geom_vline(xintercept=dat2$meany,colour="black",linetype="longdash") + 
  geom_hline(yintercept=dat2$meanWy,colour="black",linetype="longdash")+ 
  xlab("HR90 (centered on mean)") + ylab("spatial lag HR90 (centered on mean)")+
  theme(axis.line=element_line(color="black"),
        axis.title.x=element_text(size=10,vjust=0.1),
        axis.title.y=element_text(size=10,vjust=0.1),
        axis.text= element_text(colour="black", size=10, angle=0,face = "plain"),
        #plot.title = element_text(size = 10, lineheight=.8, face="bold", vjust=1), # make title bold and add space
        legend.title = element_text(size = 8), # legend title size
        legend.text = element_text(size = 8), # legend label size
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(), # removes background grid
        panel.background = element_blank(), # removes grey background; could add black axis lines with #axis.line = element_line(colour = "black")) 
        panel.spacing=unit(c(0,0,0,0), "lines"),
        plot.margin=unit(c(0,0,0,0), "mm"))  # sets margin around full plot at top, right, bottom, and left; units can also be "lines" or "cm"
dev.off()


########################################################################
#####
##### LISA analysis of DV (y) and Error (residuals)
##### this is code for Figures 2 and 3 in paper
#####
########################################################################

# DV

# Raw HR90 greyscale

plot1.1 <- ggplot(shp.df) + 
  aes(long,lat,group=group,fill=HR90) + 
  geom_polygon(colour="transparent", fill="white")
#plot1.1

plot1.1b <- plot1.1 + 
  geom_polygon(aes(long, lat, group=group), data = shp.df) +            
  #             colour = 'white', fill = 'black', alpha = .4, size = .3) +
  #scale_fill_distiller(name="LISA values", palette = "Reds", 
  #                    trans="reverse", 
  # breaks = pretty_breaks(n = 5), 
  #                   na.value="white", guide="colourbar") + 
  #scale_fill_gradient(name="LISA values", values=c("white", "red", "blue","lightblue","pink"), breaks = c("0", "1", "2", "3", "4"), labels=c("n.s.", "high-high", "low-low","low-high","hgh-low")) + 
  scale_fill_gradientn(name="HR90", 
                       #values=rescale(c(0, 20, 60)), 
                       colours=c("white","black")) +
  labs(x="", y="", title="Homicide Rates (HR90)") +
  theme(axis.ticks.y = element_blank(),axis.text.y = element_blank(), # get rid of x ticks/text
        axis.ticks.x = element_blank(),axis.text.x = element_blank(), # get rid of y ticks/text
        axis.line.x = element_blank(), axis.line.y = element_blank(), # get rid of axis lines, if present
        plot.title = element_text(size = 10, lineheight=.8, face="bold", vjust=1), # make title bold and add space
        legend.title = element_text(size = 8), # legend title size
        legend.text = element_text(size = 8), # legend label size
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(), # removes background grid
        panel.background = element_blank(), # removes grey background; could add black axis lines with #axis.line = element_line(colour = "black")) 
        panel.spacing=unit(c(0,0,0,0), "lines"),
        plot.margin=unit(c(0,0,0,0), "mm")) + # sets margin around full plot at top, right, bottom, and left; units can also be "lines" or "cm"
  coord_equal(ratio=1) +
  geom_polygon(data=shp.st, aes(long,lat, group=group), colour="grey50", fill=NA) 
#plot1.1b

# if want to generate this plot on its own, uncomment next 3 lines
#tiff(file="./figures/fig2_hr90values-bw.tif", height=6, width=9, units="in", res=300) # could also do pdf, jpeg, bmp, tiff
#print(plot1.1b)
#dev.off()

# Figure 2 LISA maps

# black and white of pos and neg lisa values with permutation method in python

plot1p <- ggplot(shp.df) + 
  aes(long,lat,group=group,fill=ls.py.I) + 
  geom_polygon(colour="transparent", fill="white")
#plot1p

plot1ppos <- plot1p + 
  geom_polygon(aes(x = long, y = lat, group = group), data = subset(shp.df, ls.py.I>0)) +            
  #             colour = 'white', fill = 'black', alpha = .4, size = .3) +
  #scale_fill_distiller(name="LISA values", palette = "Reds", 
  #                    trans="reverse", 
  # breaks = pretty_breaks(n = 5), 
  #                   na.value="white", guide="colourbar") + 
  #scale_fill_gradient(name="LISA values", values=c("white", "red", "blue","lightblue","pink"), breaks = c("0", "1", "2", "3", "4"), labels=c("n.s.", "high-high", "low-low","low-high","hgh-low")) + 
  #scale_fill_gradientn(name="LISA values", values=rescale(c(12, 4, 0, -2, -6)), colours=c("red","white","blue")) +
  scale_fill_gradientn(name="LISA values", colours=c("lightgrey","black")) +
  labs(x="", y="", title="LISA Values, positive (DV: HR90)") +
  theme(axis.ticks.y = element_blank(),axis.text.y = element_blank(), # get rid of x ticks/text
        axis.ticks.x = element_blank(),axis.text.x = element_blank(), # get rid of y ticks/text
        axis.line.x = element_blank(), axis.line.y = element_blank(), # get rid of axis lines, if present
        plot.title = element_text(size = 10, lineheight=.8, face="bold", vjust=1), # make title bold and add space
        legend.title = element_text(size = 8), # legend title size
        legend.text = element_text(size = 8), # legend label size
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(), # removes background grid
        panel.background = element_blank(), # removes grey background; could add black axis lines with #axis.line = element_line(colour = "black")) 
        panel.spacing=unit(c(0,0,0,0), "lines"),
        plot.margin=unit(c(0,0,0,0), "mm")) + # sets margin around full plot at top, right, bottom, and left; units can also be "lines" or "cm"
  coord_equal(ratio=1) +
  geom_polygon(data=shp.st, aes(long,lat, group=group), colour="grey50", fill=NA)

#plot1ppos

# if want to generate this plot on its own, uncomment next 3 lines
#tiff(file="./figures/fig2_lisavalues-perm-bw-pos.tif", height=6, width=9, units="in", res=300) # could also do pdf, jpeg, bmp, tiff
#print(plot1ppos)
#dev.off()


plot1pneg <- plot1p + 
  geom_polygon(aes(x = long, y = lat, group = group), data = subset(shp.df, ls.py.I<0)) +            
  #             colour = 'white', fill = 'black', alpha = .4, size = .3) +
  #scale_fill_distiller(name="LISA values", palette = "Reds", 
  #                    trans="reverse", 
  # breaks = pretty_breaks(n = 5), 
  #                   na.value="white", guide="colourbar") + 
  #scale_fill_gradientn(name="LISA values", values=rescale(c(12, 4, 0, -2, -6)), colours=c("red","white","blue")) +
  #scale_fill_distiller(name=NULL, palette = "Greys", direction=-1, breaks = pretty_breaks(n = 4), na.value="white", guide="colourbar") + 
  scale_fill_gradientn(name="LISA values", colours=c("black", "lightgrey")) +
  labs(x="", y="", title="LISA Values, negative (DV: HR90)") +
  theme(axis.ticks.y = element_blank(),axis.text.y = element_blank(), # get rid of x ticks/text
        axis.ticks.x = element_blank(),axis.text.x = element_blank(), # get rid of y ticks/text
        axis.line.x = element_blank(), axis.line.y = element_blank(), # get rid of axis lines, if present
        plot.title = element_text(size = 10, lineheight=.8, face="bold", vjust=1), # make title bold and add space
        legend.title = element_text(size = 8), # legend title size
        legend.text = element_text(size = 8), # legend label size
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(), # removes background grid
        panel.background = element_blank(), # removes grey background; could add black axis lines with #axis.line = element_line(colour = "black")) 
        panel.spacing=unit(c(0,0,0,0), "lines"),
        plot.margin=unit(c(0,0,0,0), "mm")) + # sets margin around full plot at top, right, bottom, and left; units can also be "lines" or "cm"
  coord_equal(ratio=1) +
  geom_polygon(data=shp.st, aes(long,lat, group=group), colour="grey50", fill=NA)

#plot1pneg

# if want to generate this plot on its own, uncomment next 3 lines
#tiff(file="./figures/fig2_lisavalues-perm-bw-neg.tif", height=6, width=9, units="in", res=300) # could also do pdf, jpeg, bmp, tiff
#print(plot1pneg)
#dev.off()

# LISA p-values with permutation method in python


# black and white p-value with permutation approach in python

plot2p <- ggplot(shp.df) + 
  aes(long,lat,group=group,fill=ls.py.pperm) + 
  geom_polygon(colour="transparent", fill="white")
#plot2p

plot2pb <- plot2p + 
  geom_polygon(aes(x = long, y = lat, group = group), data = subset(shp.df, ls.py.pperm<.2)) +            
  #             colour = 'white', fill = 'black', alpha = .4, size = .3) +
  scale_fill_gradientn(name="LISA p", 
                       #values=rescale(c(0, .05, .2)), 
                       colours=c("black", "white"), guide="colourbar") + 
  #scale_fill_gradient(name="LISA sig (p)", low="white", high="darkred", guide="colourbar") + 
  labs(x="", y="", title="LISA p-values (DV: HR90)") +
  theme(axis.ticks.y = element_blank(),axis.text.y = element_blank(), # get rid of x ticks/text
        axis.ticks.x = element_blank(),axis.text.x = element_blank(), # get rid of y ticks/text
        axis.line.x = element_blank(), axis.line.y = element_blank(), # get rid of axis lines, if present
        plot.title = element_text(size = 10, lineheight=.8, face="bold", vjust=1), # make title bold and add space
        legend.title = element_text(size = 8), # legend title size
        legend.text = element_text(size = 8), # legend label size
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(), # removes background grid
        panel.background = element_blank(), # removes grey background; could add black axis lines with #axis.line = element_line(colour = "black")) 
        panel.spacing=unit(c(0,0,0,0), "lines"),
        plot.margin=unit(c(0,0,0,0), "mm")) + # sets margin around full plot at top, right, bottom, and left; units can also be "lines" or "cm"
  coord_equal(ratio=1) +
  geom_polygon(data=shp.st, aes(long,lat, group=group), colour="grey50", fill=NA) 
#plot2pb

# if want to generate this plot on its own, uncomment next 3 lines
#tiff(file="./figures/fig2_lisap-perm-bw.tif", height=6, width=9, units="in", res=300) # could also do pdf, jpeg, bmp, tiff
#print(plot2pb)
#dev.off()


# LISA clusters with permutation approach in python - BW

plot3pyperm <- ggplot(shp.df) + 
  aes(long,lat,group=group,fill=as.factor(ls.py.clper)) + 
  geom_polygon(colour="transparent", fill="white")
#plot3perm

plot3pypermb <- plot3pyperm + 
  geom_polygon(aes(x = long, y = lat, group = group), data = shp.df) +            
  #             colour = 'white', fill = 'black', alpha = .4, size = .3) +
  #scale_fill_distiller(name=NULL, palette = "Greys", trans="reverse", breaks = pretty_breaks(n = 5), na.value="white", guide="colourbar") + 
  scale_fill_manual(name="LISA cluster", values=c("white", "black", "darkgrey","lightgrey","lightgrey"), breaks = c("0", "1", "2", "3", "4"), 
                    labels=c("n.s.", "high-high", "low-low","low-high","high-low"), guide="legend") + 
  labs(x="", y="", title="LISA Clusters (DV: HR90)") +
  theme(axis.ticks.y = element_blank(),axis.text.y = element_blank(), # get rid of x ticks/text
        axis.ticks.x = element_blank(),axis.text.x = element_blank(), # get rid of y ticks/text
        axis.line.x = element_blank(), axis.line.y = element_blank(), # get rid of axis lines, if present
        plot.title = element_text(size = 10, lineheight=.8, face="bold", vjust=1), # make title bold and add space
        legend.title = element_text(size = 8), # legend title size
        legend.text = element_text(size = 8), # legend label size
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(), # removes background grid
        panel.background = element_blank(), # removes grey background; could add black axis lines with #axis.line = element_line(colour = "black")) 
        panel.spacing=unit(c(0,0,0,0), "lines"),
        plot.margin=unit(c(0,0,0,0), "mm")) + # sets margin around full plot at top, right, bottom, and left; units can also be "lines" or "cm"
  coord_equal(ratio=1) +
  geom_polygon(data=shp.st, aes(long,lat, group=group), colour="grey50", fill=NA) 
#plot3pypermb

# if want to generate this plot on its own, uncomment next 3 lines
#tiff(file="./figures/fig2_lisa-clpermpy-bw-hh-ll.tif", height=6, width=9, units="in", res=300) # could also do pdf, jpeg, bmp, tiff
#print(plot3pypermb)
#dev.off()

plot3pypermc <- plot3pyperm + 
  geom_polygon(aes(x = long, y = lat, group = group), data = shp.df) +            
  #             colour = 'white', fill = 'black', alpha = .4, size = .3) +
  #scale_fill_distiller(name=NULL, palette = "Greys", trans="reverse", breaks = pretty_breaks(n = 5), na.value="white", guide="colourbar") + 
  scale_fill_manual(name="LISA cluster", values=c("white", "lightgrey", "lightgrey","darkgrey","black"), breaks = c("0", "1", "2", "3", "4"), 
                    labels=c("n.s.", "high-high", "low-low","low-high","high-low"), guide="legend") + 
  labs(x="", y="", title="LISA Clusters (DV: HR90)") +
  theme(axis.ticks.y = element_blank(),axis.text.y = element_blank(), # get rid of x ticks/text
        axis.ticks.x = element_blank(),axis.text.x = element_blank(), # get rid of y ticks/text
        axis.line.x = element_blank(), axis.line.y = element_blank(), # get rid of axis lines, if present
        plot.title = element_text(size = 10, lineheight=.8, face="bold", vjust=1), # make title bold and add space
        legend.title = element_text(size = 8), # legend title size
        legend.text = element_text(size = 8), # legend label size
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(), # removes background grid
        panel.background = element_blank(), # removes grey background; could add black axis lines with #axis.line = element_line(colour = "black")) 
        panel.spacing=unit(c(0,0,0,0), "lines"),
        plot.margin=unit(c(0,0,0,0), "mm")) + # sets margin around full plot at top, right, bottom, and left; units can also be "lines" or "cm"
  coord_equal(ratio=1) +
  geom_polygon(data=shp.st, aes(long,lat, group=group), colour="grey50", fill=NA) 
#plot3pypermc

# if want to generate this plot on its own, uncomment next 3 lines
#tiff(file="./figures/fig2_lisa-clpermpy-bw-lh-hl.tif", height=6, width=9, units="in", res=300) # could also do pdf, jpeg, bmp, tiff
#print(plot3pypermc)
#dev.off()

###################################################
# 
# Figure 2: Combined Graphs with Permutation Approach - BW
#
###################################################

tiff(file="./figures/fig2_comb_pyperm_bw.tif", height=6, width=9, units="in", res=300) # could also do pdf, jpeg, bmp, tiff
grid.arrange(plot1.1b, plot2pb, plot1ppos, plot1pneg, plot3pypermb, plot3pypermc, ncol=2)
dev.off()


#######################################
#
# Code for Figure 3 with Residuals 
# Repeat maps with residual LISA stats
#
#######################################

# Residuals: ols2res

plot4 <- ggplot(shp.df) + 
  aes(long,lat,group=group,fill=ols2res) + 
  geom_polygon(colour="transparent", fill="white")
#plot4

# black and white of pos and neg residuals

plot4a.pos <- plot4 + 
  geom_polygon(aes(x = long, y = lat, group = group), data = subset(shp.df, ols2res>0)) +            
  #             colour = 'white', fill = 'black', alpha = .4, size = .3) +
  #scale_fill_distiller(name="LISA values", palette = "Reds", 
  #                    trans="reverse", 
  # breaks = pretty_breaks(n = 5), 
  #                   na.value="white", guide="colourbar") + 
  #scale_fill_gradient(name="LISA values", values=c("white", "red", "blue","lightblue","pink"), breaks = c("0", "1", "2", "3", "4"), labels=c("n.s.", "high-high", "low-low","low-high","hgh-low")) + 
  #scale_fill_gradientn(name="LISA values", values=rescale(c(12, 4, 0, -2, -6)), colours=c("red","white","blue")) +
  scale_fill_gradientn(name="Residuals", colours=c("lightgrey","black")) +
  labs(x="", y="", title="OLS Residuals (positive)") +
  theme(axis.ticks.y = element_blank(),axis.text.y = element_blank(), # get rid of x ticks/text
        axis.ticks.x = element_blank(),axis.text.x = element_blank(), # get rid of y ticks/text
        axis.line.x = element_blank(), axis.line.y = element_blank(), # get rid of axis lines, if present
        plot.title = element_text(size = 10, lineheight=.8, face="bold", vjust=1), # make title bold and add space
        legend.title = element_text(size = 8), # legend title size
        legend.text = element_text(size = 8), # legend label size
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(), # removes background grid
        panel.background = element_blank(), # removes grey background; could add black axis lines with #axis.line = element_line(colour = "black")) 
        panel.spacing=unit(c(0,0,0,0), "lines"),
        plot.margin=unit(c(0,0,0,0), "mm")) + # sets margin around full plot at top, right, bottom, and left; units can also be "lines" or "cm"
  coord_equal(ratio=1) +
  geom_polygon(data=shp.st, aes(long,lat, group=group), colour="grey50", fill=NA)
#plot4a.pos

# if want to generate this plot on its own, uncomment next 3 lines
#tiff(file="./figures/fig3_olsresiduals-bw-pos.tif", height=6, width=9, units="in", res=300) # could also do pdf, jpeg, bmp, tiff
#print(plot4a.pos)
#dev.off()


plot4a.neg <- plot4 + 
  geom_polygon(aes(x = long, y = lat, group = group), data = subset(shp.df, ols2res<0)) +            
  #             colour = 'white', fill = 'black', alpha = .4, size = .3) +
  #scale_fill_distiller(name="LISA values", palette = "Reds", 
  #                    trans="reverse", 
  # breaks = pretty_breaks(n = 5), 
  #                   na.value="white", guide="colourbar") + 
  #scale_fill_gradient(name="LISA values", values=c("white", "red", "blue","lightblue","pink"), breaks = c("0", "1", "2", "3", "4"), labels=c("n.s.", "high-high", "low-low","low-high","hgh-low")) + 
  #scale_fill_gradientn(name="LISA values", values=rescale(c(12, 4, 0, -2, -6)), colours=c("red","white","blue")) +
  scale_fill_gradientn(name="Residuals", colours=c("black","lightgrey")) +
  labs(x="", y="", title="OLS Residuals (negative)") +
  theme(axis.ticks.y = element_blank(),axis.text.y = element_blank(), # get rid of x ticks/text
        axis.ticks.x = element_blank(),axis.text.x = element_blank(), # get rid of y ticks/text
        axis.line.x = element_blank(), axis.line.y = element_blank(), # get rid of axis lines, if present
        plot.title = element_text(size = 10, lineheight=.8, face="bold", vjust=1), # make title bold and add space
        legend.title = element_text(size = 8), # legend title size
        legend.text = element_text(size = 8), # legend label size
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(), # removes background grid
        panel.background = element_blank(), # removes grey background; could add black axis lines with #axis.line = element_line(colour = "black")) 
        panel.spacing=unit(c(0,0,0,0), "lines"),
        plot.margin=unit(c(0,0,0,0), "mm")) + # sets margin around full plot at top, right, bottom, and left; units can also be "lines" or "cm"
  coord_equal(ratio=1) +
  geom_polygon(data=shp.st, aes(long,lat, group=group), colour="grey50", fill=NA)
#plot4a.neg

# if want to generate this plot on its own, uncomment next 3 lines
#tiff(file="./figures/fig3_olsresiduals-bw-neg.tif", height=6, width=9, units="in", res=300) # could also do pdf, jpeg, bmp, tiff
#print(plot4a.neg)
#dev.off()


#######################################################
# LISA of residuals using permutation method in Python
#######################################################

#names(shp.df)

# raw residual plots above
# color:
# bw: 

# lisa values of residuals using permutation approach in python - BW

plot4.2pyperm <- ggplot(shp.df) + 
  aes(long,lat,group=group,fill=ls.py.lresp) + 
  geom_polygon(colour="transparent", fill="white")
#plot4.2pyperm

plot4bpyperm.pos <- plot4.2pyperm + 
  geom_polygon(aes(x = long, y = lat, group = group), data = subset(shp.df, ls.py.lresp>0)) +            
  #             colour = 'white', fill = 'black', alpha = .4, size = .3) +
  #scale_fill_distiller(name="LISA values", palette = "Reds", 
  #                    trans="reverse", 
  # breaks = pretty_breaks(n = 5), 
  #                   na.value="white", guide="colourbar") + 
  #scale_fill_gradient(name="LISA values", values=c("white", "red", "blue","lightblue","pink"), breaks = c("0", "1", "2", "3", "4"), labels=c("n.s.", "high-high", "low-low","low-high","hgh-low")) + 
  #scale_fill_gradientn(name="LISA values", values=rescale(c(12, 4, 0, -2, -6)), colours=c("red","white","blue")) +
  scale_fill_gradientn(name="LISA values", colours=c("lightgrey","black")) +
  labs(x="", y="", title="LISA Values, positive (residuals)") +
  theme(axis.ticks.y = element_blank(),axis.text.y = element_blank(), # get rid of x ticks/text
        axis.ticks.x = element_blank(),axis.text.x = element_blank(), # get rid of y ticks/text
        axis.line.x = element_blank(), axis.line.y = element_blank(), # get rid of axis lines, if present
        plot.title = element_text(size = 10, lineheight=.8, face="bold", vjust=1), # make title bold and add space
        legend.title = element_text(size = 8), # legend title size
        legend.text = element_text(size = 8), # legend label size
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(), # removes background grid
        panel.background = element_blank(), # removes grey background; could add black axis lines with #axis.line = element_line(colour = "black")) 
        panel.spacing=unit(c(0,0,0,0), "lines"),
        plot.margin=unit(c(0,0,0,0), "mm")) + # sets margin around full plot at top, right, bottom, and left; units can also be "lines" or "cm"
  coord_equal(ratio=1) +
  geom_polygon(data=shp.st, aes(long,lat, group=group), colour="grey50", fill=NA)
#plot4bpyperm.pos

# if want to generate this plot on its own, uncomment next 3 lines
#tiff(file="./figures/fig3_lisavalues-pyperm-res-bw-pos.tif", height=6, width=9, units="in", res=300) # could also do pdf, jpeg, bmp, tiff
#print(plot4bpyperm.pos)
#dev.off()

plot4bpyperm.neg <- plot4.2pyperm + 
  geom_polygon(aes(x = long, y = lat, group = group), data = subset(shp.df, ls.py.lresp<0)) +            
  #             colour = 'white', fill = 'black', alpha = .4, size = .3) +
  #scale_fill_distiller(name="LISA values", palette = "Reds", 
  #                    trans="reverse", 
  # breaks = pretty_breaks(n = 5), 
  #                   na.value="white", guide="colourbar") + 
  #scale_fill_gradient(name="LISA values", values=c("white", "red", "blue","lightblue","pink"), breaks = c("0", "1", "2", "3", "4"), labels=c("n.s.", "high-high", "low-low","low-high","hgh-low")) + 
  #scale_fill_gradientn(name="LISA values", values=rescale(c(12, 4, 0, -2, -6)), colours=c("red","white","blue")) +
  scale_fill_gradientn(name="LISA values", colours=c("black","lightgrey")) +
  labs(x="", y="", title="LISA Values, negative (residuals)") +
  theme(axis.ticks.y = element_blank(),axis.text.y = element_blank(), # get rid of x ticks/text
        axis.ticks.x = element_blank(),axis.text.x = element_blank(), # get rid of y ticks/text
        axis.line.x = element_blank(), axis.line.y = element_blank(), # get rid of axis lines, if present
        plot.title = element_text(size = 10, lineheight=.8, face="bold", vjust=1), # make title bold and add space
        legend.title = element_text(size = 8), # legend title size
        legend.text = element_text(size = 8), # legend label size
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(), # removes background grid
        panel.background = element_blank(), # removes grey background; could add black axis lines with #axis.line = element_line(colour = "black")) 
        panel.spacing=unit(c(0,0,0,0), "lines"),
        plot.margin=unit(c(0,0,0,0), "mm")) + # sets margin around full plot at top, right, bottom, and left; units can also be "lines" or "cm"
  coord_equal(ratio=1) +
  geom_polygon(data=shp.st, aes(long,lat, group=group), colour="grey50", fill=NA)
#plot4bpyperm.neg

# if want to generate this plot on its own, uncomment next 3 lines
#tiff(file="./figures/fig3_lisavalues-pyperm-res-bw-neg.tif", height=6, width=9, units="in", res=300) # could also do pdf, jpeg, bmp, tiff
#print(plot4bpyperm.neg)
#dev.off()

# black and white permutation approach in python

plot5pyperm <- ggplot(shp.df) + 
  aes(long,lat,group=group,fill=ls.py.presp) + 
  geom_polygon(colour="transparent", fill="white")
#plot5pyperm

plot5pypermb <- plot5pyperm + 
  geom_polygon(aes(x = long, y = lat, group = group), data = subset(shp.df, ls.py.presp<.2)) +            
  #             colour = 'white', fill = 'black', alpha = .4, size = .3) +
  scale_fill_gradientn(name="LISA p", 
                       #values=rescale(c(0, .05, .2)), 
                       colours=c("black", "white"), guide="colourbar") + 
  #scale_fill_gradient(name="LISA sig (p)", low="white", high="darkred", guide="colourbar") + 
  labs(x="", y="", title="LISA p-values (residuals)") +
  theme(axis.ticks.y = element_blank(),axis.text.y = element_blank(), # get rid of x ticks/text
        axis.ticks.x = element_blank(),axis.text.x = element_blank(), # get rid of y ticks/text
        axis.line.x = element_blank(), axis.line.y = element_blank(), # get rid of axis lines, if present
        plot.title = element_text(size = 10, lineheight=.8, face="bold", vjust=1), # make title bold and add space
        legend.title = element_text(size = 8), # legend title size
        legend.text = element_text(size = 8), # legend label size
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(), # removes background grid
        panel.background = element_blank(), # removes grey background; could add black axis lines with #axis.line = element_line(colour = "black")) 
        panel.spacing=unit(c(0,0,0,0), "lines"),
        plot.margin=unit(c(0,0,0,0), "mm")) + # sets margin around full plot at top, right, bottom, and left; units can also be "lines" or "cm"
  coord_equal(ratio=1) +
  geom_polygon(data=shp.st, aes(long,lat, group=group), colour="grey50", fill=NA) 
#plot5pypermb

# if want to generate this plot on its own, uncomment next 3 lines
#tiff(file="./figures/fig3_lisap-pyperm-res-bw.tif", height=6, width=9, units="in", res=300) # could also do pdf, jpeg, bmp, tiff
#print(plot5pypermb)
#dev.off()

# LISA clusters - residuals using permutation method in python - BW

plot6pyperm <- ggplot(shp.df) + 
  aes(long,lat,group=group,fill=as.factor(ls.py.qresp)) + 
  geom_polygon(colour="transparent", fill="white")
#plot6pyperm

plot6pypermb <- plot6pyperm + 
  geom_polygon(aes(x = long, y = lat, group = group), data = shp.df) +            
  #             colour = 'white', fill = 'black', alpha = .4, size = .3) +
  #scale_fill_distiller(name=NULL, palette = "Greys", trans="reverse", breaks = pretty_breaks(n = 5), na.value="white", guide="colourbar") + 
  scale_fill_manual(name="LISA cluster", values=c("white", "black", "darkgrey","lightgrey","lightgrey"), breaks = c("0", "1", "2", "3", "4"), 
                    labels=c("n.s.", "high-high", "low-low","low-high","high-low"), guide="legend") + 
  labs(x="", y="", title="LISA Clusters (residuals)") +
  theme(axis.ticks.y = element_blank(),axis.text.y = element_blank(), # get rid of x ticks/text
        axis.ticks.x = element_blank(),axis.text.x = element_blank(), # get rid of y ticks/text
        axis.line.x = element_blank(), axis.line.y = element_blank(), # get rid of axis lines, if present
        plot.title = element_text(size = 10, lineheight=.8, face="bold", vjust=1), # make title bold and add space
        legend.title = element_text(size = 8), # legend title size
        legend.text = element_text(size = 8), # legend label size
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(), # removes background grid
        panel.background = element_blank(), # removes grey background; could add black axis lines with #axis.line = element_line(colour = "black")) 
        panel.spacing=unit(c(0,0,0,0), "lines"),
        plot.margin=unit(c(0,0,0,0), "mm")) + # sets margin around full plot at top, right, bottom, and left; units can also be "lines" or "cm"
  coord_equal(ratio=1) +
  geom_polygon(data=shp.st, aes(long,lat, group=group), colour="grey50", fill=NA) 
#plot6pypermb

# if want to generate this plot on its own, uncomment next 3 lines
#tiff(file="./figures/fig3_lisacl-permpy-res-bw-hh-ll.tif", height=6, width=9, units="in", res=300) # could also do pdf, jpeg, bmp, tiff
#print(plot6pypermb)
#dev.off()

plot6pypermc <- plot6pyperm + 
  geom_polygon(aes(x = long, y = lat, group = group), data = shp.df) +            
  #             colour = 'white', fill = 'black', alpha = .4, size = .3) +
  #scale_fill_distiller(name=NULL, palette = "Greys", trans="reverse", breaks = pretty_breaks(n = 5), na.value="white", guide="colourbar") + 
  scale_fill_manual(name="LISA cluster", values=c("white", "lightgrey", "lightgrey","darkgrey","black"), breaks = c("0", "1", "2", "3", "4"), 
                    labels=c("n.s.", "high-high", "low-low","low-high","high-low"), guide="legend") + 
  labs(x="", y="", title="LISA Clusters (residuals)") +
  theme(axis.ticks.y = element_blank(),axis.text.y = element_blank(), # get rid of x ticks/text
        axis.ticks.x = element_blank(),axis.text.x = element_blank(), # get rid of y ticks/text
        axis.line.x = element_blank(), axis.line.y = element_blank(), # get rid of axis lines, if present
        plot.title = element_text(size = 10, lineheight=.8, face="bold", vjust=1), # make title bold and add space
        legend.title = element_text(size = 8), # legend title size
        legend.text = element_text(size = 8), # legend label size
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(), # removes background grid
        panel.background = element_blank(), # removes grey background; could add black axis lines with #axis.line = element_line(colour = "black")) 
        panel.spacing=unit(c(0,0,0,0), "lines"),
        plot.margin=unit(c(0,0,0,0), "mm")) + # sets margin around full plot at top, right, bottom, and left; units can also be "lines" or "cm"
  coord_equal(ratio=1) +
  geom_polygon(data=shp.st, aes(long,lat, group=group), colour="grey50", fill=NA) 
#plot6pypermc

# if want to generate this plot on its own, uncomment next 3 lines
#tiff(file="./figures/fig3_lisacl-permpy-res-bw-lh-hl.tif", height=6, width=9, units="in", res=300) # could also do pdf, jpeg, bmp, tiff
#print(plot6pypermc)
#dev.off()


###################################################
# Figure 3: Combined Graphs with Permutation Approach - BW
###################################################

tiff(file="./figures/fig3_comb_pyperm_bw.tif", height=6, width=9, units="in", res=300) # could also do pdf, jpeg, bmp, tiff
grid.arrange(plot4a.pos, plot4a.neg, plot4bpyperm.pos, plot4bpyperm.neg, 
             plot5pypermb, nullGrob(), plot6pypermb, plot6pypermc, ncol=2)
dev.off()


###############################
#
# TABLE 2
#
###############################

####################################################
#
# Generate spatial lag of DV
#
####################################################

shp$HR90lag <- lag.listw(w, shp$HR90)

# high-high

tab2 <- subset(shp@data, 
               ls.py.clper==1,
               select=c(FIPSNO, NAME, STATE_NAME, HR90, HR90lag, ls.py.I, ls.py.pperm))
names(tab2)

tab2cl1 <- tab2[order(-tab2$ls.py.I, tab2$ls.py.pperm), ]
tab2cl1[1:15,]
str(tab2cl1)

tab2cl1b <- data.frame(FIPSNO=integer(), NAME=character(), STATE_NAME=character(), 
                       HR90=numeric(), HR90lag=numeric(), 
                       ls.py.I=numeric(), ls.py.pperm=numeric())
str(tab2cl1b)

shpdata <- shp@data
shpdata <- subset(shpdata, ls.py.clper==1)
str(shpdata)
# should show an object with 739 rows and 80 cols

for (i in 1:length(shpdata[,1])) {
  if (tab2cl1[i,7] < .05) {
    tab2cl1b <- rbind(tab2cl1b, tab2cl1[i,])
  }
}

tab2cl1b[1:5,]

# should show:
#FIPSNO         NAME     STATE_NAME     HR90  HR90lag ls.py.I ls.py.pperm
#28055  28055    Issaquena Mississippi 34.92230 23.05810  10.397       0.002
#22035  22035 East Carroll   Louisiana 37.76565 19.21942   9.335       0.002
#28151  28151   Washington Mississippi 30.42124 23.12162   9.286       0.002
#28125  28125      Sharkey Mississippi 28.30456 22.80082   8.334       0.002
#28083  28083      Leflore Mississippi 27.67289 20.62248   6.975       0.002


# low-low clusters
tab2 <- subset(shp@data, 
               ls.py.clper==2,
               select=c(FIPSNO, NAME, STATE_NAME, HR90, HR90lag, ls.py.I, ls.py.pperm))

tab2cl2 <- tab2[order(-tab2$ls.py.I, tab2$ls.py.pperm), ]
tab2cl2[1:15,]
str(tab2cl2)

tab2cl2b <- data.frame(FIPSNO=integer(), NAME=character(), STATE_NAME=character(), 
                       HR90=numeric(), HR90lag=numeric(), 
                       ls.py.I=numeric(), ls.py.pperm=numeric())
str(tab2cl2b)

shpdata <- shp@data
shpdata <- subset(shpdata, ls.py.clper==2)
str(shpdata)
# should show an object with 739 rows and 80 cols

for (i in 1:length(shpdata[,1])) {
  if (tab2cl2[i,7] < .05) {
    tab2cl2b <- rbind(tab2cl2b, tab2cl2[i,])
  }
}

tab2cl2b[1:5,]

# should show:
#      FIPSNO      NAME   STATE_NAME HR90 HR90lag ls.py.I ls.py.pperm
#705   31101     Keith     Nebraska    0  0.0854   0.867       0.002
#787   31135   Perkins     Nebraska    0  0.0789   0.867       0.002
#1211  20109     Logan       Kansas    0  0.0000   0.867       0.002
#1618  20007    Barber       Kansas    0  0.3194   0.867       0.002
#202   46089 McPherson South Dakota    0  0.0852   0.854       0.002

# low-high

tab2 <- subset(shp@data, 
               ls.py.clper==3,
               select=c(FIPSNO, NAME, STATE_NAME, HR90, HR90lag, ls.py.I, ls.py.pperm))

tab2cl3 <- tab2[order(tab2$ls.py.I, tab2$ls.py.pperm), ]
tab2cl3[1:15,]
str(tab2cl3)

tab2cl3b <- data.frame(FIPSNO=integer(), NAME=character(), STATE_NAME=character(), 
                       HR90=numeric(), HR90lag=numeric(), 
                       ls.py.I=numeric(), ls.py.pperm=numeric())
str(tab2cl3b)

shpdata <- shp@data
shpdata <- subset(shpdata, ls.py.clper==3)
str(shpdata)
# should show an object with 84 rows and 80 cols

for (i in 1:length(shpdata[,1])) {
  if (tab2cl3[i,7] < .05) {
    tab2cl3b <- rbind(tab2cl3b, tab2cl3[i,])
  }
}

tab2cl3b[1:5,]

# should show
#FIPSNO      NAME  STATE_NAME HR90  HR90lag ls.py.I ls.py.pperm
#13183  13183      Long  Georgia    0 13.08078  -1.054       0.010
#8007    8007 Archuleta Colorado    0 11.62266  -1.013       0.008
#48365  48365    Panola    Texas    0 12.65118  -0.947       0.018
#21109  21109   Jackson Kentucky    0 12.72195  -0.917       0.012
#48271  48271    Kinney    Texas    0 12.62916  -0.904       0.016

# high-low

tab2 <- subset(shp@data, 
               ls.py.clper==4,
               select=c(FIPSNO, NAME, STATE_NAME, HR90, HR90lag, ls.py.I, ls.py.pperm))

tab2cl4 <- tab2[order(tab2$ls.py.I, tab2$ls.py.pperm), ]
tab2cl4[1:15,]
str(tab2cl4)

tab2cl4b <- data.frame(FIPSNO=integer(), NAME=character(), STATE_NAME=character(), 
                       HR90=numeric(), HR90lag=numeric(), 
                       ls.py.I=numeric(), ls.py.pperm=numeric())
str(tab2cl4b)

shpdata <- shp@data
shpdata <- subset(shpdata, ls.py.clper==4)
str(shpdata)
# should show an object with 70 rows and 80 cols

for (i in 1:length(shpdata[,1])) {
  if (tab2cl4[i,7] < .05) {
    tab2cl4b <- rbind(tab2cl4b, tab2cl4[i,])
  }
}

tab2cl4b[1:5,]

# should show
#FIPSNO       NAME      STATE_NAME     HR90   HR90lag ls.py.I ls.py.pperm
#8053    8053   Hinsdale     Colorado 71.37759 3.5810709  -6.830       0.004
#46063  46063    Harding South Dakota 19.97204 0.2934208  -1.832       0.002
#31173  31173   Thurston     Nebraska 19.22338 1.1480943  -1.455       0.004
#17035  17035 Cumberland     Illinois 18.74414 1.8367758  -1.282       0.002
#41023  41023      Grant       Oregon 16.97865 0.9924572  -1.253       0.002

# not significant

tab2 <- subset(shp@data, 
               ls.py.clper==0,
               select=c(FIPSNO, NAME, STATE_NAME, HR90, HR90lag, ls.py.I, ls.py.pperm))

tab2cl5 <- tab2[order(-tab2$ls.py.pperm), ]
tab2cl5[1:5,]
str(tab2cl5)

# table 2 as appears in paper
table2 <- rbind(tab2cl1b[1:5,], tab2cl2b[1:5,], 
                tab2cl3b[1:5,], tab2cl4b[1:5,])

# tabl2 2 from paper plus non-significant locations
#table2 <- rbind(tab2cl1b[1:5,], tab2cl2b[1:5,], 
#                tab2cl3b[1:5,], tab2cl4b[1:5,], 
#                tab2cl5[1:5,])
table2

# write table to file
write.csv(table2, file = "./tables/table2.csv")


###################################################
#
# TABLE 3
#
###################################################

####################################################
#
# Generate spatial lag of residuals
#
####################################################

shp$ols2reslag <- lag.listw(w, shp$ols2res)

# high-high

tab3 <- subset(shp@data, 
               ls.py.qresp==1,
               select=c(FIPSNO, NAME, STATE_NAME, 
                        ols2res, ols2reslag, 
                        ls.py.lresp, ls.py.presp))
names(tab3)

tab3cl1 <- tab3[order(-tab3$ls.py.lresp, tab3$ls.py.presp), ]
cbind(tab3cl1[1:15,1:3], round(tab3cl1[1:15,4:7], 3))
str(tab3cl1)

tab3cl1b <- data.frame(FIPSNO=integer(), NAME=character(), STATE_NAME=character(), 
                       ols2res=numeric(), ols2reslag=numeric(), 
                       ls.py.I=numeric(), ls.py.pperm=numeric())
str(tab3cl1b)

shpdata <- shp@data
shpdata <- subset(shpdata, ls.py.qresp==1)
str(shpdata)
# should show an object with 168 rows and 80 cols


for (i in 1:length(shpdata[,1])) {
  if (tab3cl1[i,7] < .05) {
    tab3cl1b <- rbind(tab3cl1b, tab3cl1[i,])
  }
}

tab3cl1b[1:5,]

# should show
#     FIPSNO           NAME           STATE_NAME ols2res ols2reslag ls.py.lresp
#1251  11001     Washington District of Columbia    43.5       4.73        6.80
#1366  49031          Piute                 Utah    27.7       3.02        4.33
#1153  24510 Baltimore City             Maryland    13.9       6.43        3.72
#1215  24033 Prince Georges             Maryland    12.4       7.18        3.65
#2731  48435         Sutton                Texas    20.9       2.70        3.50
#          ls.py.presp
#1251       0.032
#1366       0.026
#1153       0.002
#1215       0.002
#2731       0.040

# low-low clusters
tab3 <- subset(shp@data, 
               ls.py.qresp==2,
               select=c(FIPSNO, NAME, STATE_NAME, 
                        ols2res, ols2reslag, 
                        ls.py.lresp, ls.py.presp))

tab3cl2 <- tab3[order(-tab3$ls.py.lresp, tab3$ls.py.presp), ]
cbind(tab3cl2[1:15,1:3], round(tab3cl2[1:15,4:7], 3))
str(tab3cl2)

tab3cl2b <- data.frame(FIPSNO=integer(), NAME=character(), STATE_NAME=character(), 
                       ols2res=numeric(), ols2reslag=numeric(), 
                       ls.py.I=numeric(), ls.py.pperm=numeric())
str(tab3cl2b)

shpdata <- shp@data
shpdata <- subset(shpdata, ls.py.qresp==2)
str(shpdata)
# should show an object with 170 rows and 80 cols

for (i in 1:length(shpdata[,1])) {
  if (tab3cl2[i,7] < .05) {
    tab3cl2b <- rbind(tab3cl2b, tab3cl2[i,])
  }
}

tab3cl2b[1:5,]

# should show:
#     FIPSNO      NAME  STATE_NAME ols2res ols2reslag ls.py.lresp ls.py.presp
#2455   1105     Perry     Alabama  -19.10      -4.25        2.84       0.010
#2661  28157 Wilkinson Mississippi  -12.29      -5.34        2.65       0.002
#2422   1065      Hale     Alabama  -15.10      -4.18        2.46       0.002
#2669  22009 Avoyelles   Louisiana   -8.66      -4.97        1.90       0.002
#2611  22029 Concordia   Louisiana   -8.99      -5.42        1.90       0.004

# low-high

tab3 <- subset(shp@data, 
               ls.py.qresp==3,
               select=c(FIPSNO, NAME, STATE_NAME, 
                        ols2res, ols2reslag, 
                        ls.py.lresp, ls.py.presp))

tab3cl3 <- tab3[order(-tab3$ls.py.lresp, tab3$ls.py.presp), ]
cbind(tab3cl3[1:15,1:3], round(tab3cl3[1:15,4:7], 3))
str(tab3cl3)

tab3cl3b <- data.frame(FIPSNO=integer(), NAME=character(), STATE_NAME=character(), 
                       ols2res=numeric(), ols2reslag=numeric(), 
                       ls.py.I=numeric(), ls.py.pperm=numeric())
str(tab3cl3b)

shpdata <- shp@data
shpdata <- subset(shpdata, ls.py.qresp==3)
str(shpdata)
# should show an object with 79 rows and 80 cols

for (i in 1:length(shpdata[,1])) {
  if (tab3cl3[i,7] < .05) {
    tab3cl3b <- rbind(tab3cl3b, tab3cl3[i,])
  }
}

tab3cl3b[1:5,]

# should show

#     FIPSNO      NAME     STATE_NAME ols2res ols2reslag ls.py.lresp ls.py.presp
#2821  48039  Brazoria          Texas -0.0852       3.71      -0.012       0.048
#1841  47159     Smith      Tennessee -0.1351       3.26      -0.021       0.026
#2200  45047 Greenwood South Carolina -0.1326       4.12      -0.022       0.018
#1773  47111     Macon      Tennessee -0.2947       3.45      -0.046       0.048
#2134  28137      Tate    Mississippi -0.4350       3.71      -0.063       0.042

# high-low

tab3 <- subset(shp@data, 
               ls.py.qresp==4,
               select=c(FIPSNO, NAME, STATE_NAME, 
                        ols2res, ols2reslag, 
                        ls.py.lresp, ls.py.presp))

tab3cl4 <- tab3[order(-tab3$ls.py.lresp, tab3$ls.py.presp), ]
cbind(tab3cl4[1:15,1:3], round(tab3cl4[1:15,4:7], 3))
str(tab3cl4)

tab3cl4b <- data.frame(FIPSNO=integer(), NAME=character(), STATE_NAME=character(), 
                       ols2res=numeric(), ols2reslag=numeric(), 
                       ls.py.I=numeric(), ls.py.pperm=numeric())
str(tab3cl4b)

shpdata <- shp@data
shpdata <- subset(shpdata, ls.py.qresp==4)
str(shpdata)
# should show an object with 58 rows and 80 cols

for (i in 1:length(shpdata[,1])) {
  if (tab3cl4[i,7] < .05) {
    tab3cl4b <- rbind(tab3cl4b, tab3cl4[i,])
  }
}

tab3cl4b[1:5,]

# should show

#     FIPSNO    NAME   STATE_NAME ols2res ols2reslag ls.py.lresp ls.py.presp
#2955  26051 Gladwin     Michigan  0.0214      -2.57      -0.003       0.024
#723   39103  Medina         Ohio  0.1014      -2.62      -0.011       0.034
#1471  29217  Vernon     Missouri  0.1260      -3.08      -0.015       0.020
#416   46123   Tripp South Dakota  0.1820      -3.54      -0.027       0.008
#1570  21137 Lincoln     Kentucky  0.5985      -3.27      -0.076       0.016


# not significant

tab3 <- subset(shp@data, 
               ls.py.qresp==0,
               select=c(FIPSNO, NAME, STATE_NAME, 
                        ols2res, ols2reslag, 
                        ls.py.lresp, ls.py.presp))

tab3cl5 <- tab3[order(-tab3$ls.py.presp), ]
cbind(tab3cl5[1:15,1:3], round(tab3cl5[1:15,4:7], 3))
str(tab3cl5)

# table3 as appears in paper
table3 <- rbind(tab3cl1b[1:5,], tab3cl2b[1:5,])

# table 3 from paper plus low-high, high-low, and non-significant locations
#table3 <- rbind(tab3cl1b[1:5,], tab3cl2b[1:5,], 
#                tab3cl3b[1:5,], tab3cl4b[1:5,], 
#                tab3cl5[1:5,])
table3

# write table to file
write.csv(table3, file = "./tables/table3.csv")


# save workspace
save.image("./data/working/harbersingram_psrm_20181206.RData")

##############################
#
# APPENDIX
# NOTE:
# If want to reproduce figures in appendix,
# then uncomment line below (this will run separate script
# for appendix; if want to run that script on its own, read
# header at top of that script) 
#
##############################

source('./code/harbersingram_psrm4_20181206_appendix.R')

#end